function [adj,refi,a1,b1,c] = ...
    getimpact(globvar,par,na1,rb,rb1,V_bf,V_af,b0,a0,z0,p0,f0,fmat,fprob,t) 

%t=1;n=2;i=4;rb=rb_bf;rb1=rb_af;b0=BINI(i);a0=AINI(i);z0=ZINI(i,1,2);p0=PINI(n);f0=fini(i);

s=par{1};beta=par{2};theta=par{3};ra=par{4};fb=par{5};mar=par{6};infi=par{7};fH=par{10};x=par{11};rb_af=par{12};alpha=par{13};

amat  = globvar{1};
zmat  = globvar{2}; 
ztran = globvar{3};
pmat  = globvar{4};   
ptran = globvar{5};
bmat  = (1-theta)*exp(pmat);

nz = length(zmat);
np = length(pmat);
nf = length(fmat);

% choice set of a1
a1choice  = linspace(0,max(amat),na1);

% points to be interpolated
a1interp  = repmat(a1choice',1,nz*np*nf);
b1ninterp = repmat(b0,na1,nz*np*nf);   % norefi: b'=b (same debt level)

b1r = (1-theta)*exp(p0); % refi: b'=(1-theta)p (at collateral constraint)
b1rinterp = repmat(b1r,na1,nz*np*nf);   
z1interp  = repmat(zmat',na1,np*nf);
p1interp  = repmat(repelem(pmat',1,nz),na1,nf);
f1interp  = repelem(fmat',na1,nz*np);


% transition prob
zprob = ztran(find(zmat==z0),:);
pprob = ptran(find(pmat==p0),:);
zpf_prob = kron(fprob(1,:),kron(pprob,zprob));

% consumption and u
wlthn = exp(z0)+(1+ra)*a0-(1+rb)*b0;  
wlthr = exp(z0)+(1+ra)*a0-(1+rb1)*b0;  

cr = repmat(wlthr+(1-fb)*b1r,1,na1)-a1choice;  %ns*na1
ur = u(cr,s,infi,mar) + f0 + (f0==0|f0==fH)*(rb-rb1)*100*x*(1/t^alpha);

cn = repmat(wlthn+b0,1,na1)-a1choice;    %ns*na1
un = u(cn,s,infi,mar);

% ============== value function  =============
    % Interpolate
Vn = V_bf*(rb~=rb_af)+V_af*(rb==rb_af);

if fprob==1   % frictionless case
v1nvec = interpn(bmat,amat,zmat,pmat,Vn,b1ninterp,a1interp,z1interp,p1interp,'linear');
v1rvec = interpn(bmat,amat,zmat,pmat,V_af,b1rinterp,a1interp,z1interp,p1interp,'linear');
else
v1nvec = interpn(bmat,amat,zmat,pmat,fmat,Vn,b1ninterp,a1interp,z1interp,p1interp,f1interp,'linear');
v1rvec = interpn(bmat,amat,zmat,pmat,fmat,V_af,b1rinterp,a1interp,z1interp,p1interp,f1interp,'linear');
end
v1nvec(isnan(v1nvec)) = -infi;
v1rvec(isnan(v1rvec)) = -infi;

    % expectation
ev1nvec = (v1nvec*zpf_prob')';
ev1rvec = (v1rvec*zpf_prob')';

    % value function of all states, all choices
vn  = un +beta*ev1nvec;
vr  = ur +beta*ev1rvec;

    % max Vn, max Vr
[maxvr,maxvridx] = max(vr,[],2);
[maxvn,maxvnidx] = max(vn,[],2);

adj = (maxvr>maxvn);
a1 = adj*a1choice(maxvridx)+(1-adj)*a1choice(maxvnidx);
b1 = adj*b1r+(1-adj)*b0;
c  = adj*cr(maxvridx)+(1-adj)*cn(maxvnidx);
refi = adj*(b1>b0);



